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A multivariate generalization of 
Prony's method 

Stefan Kunis^^ Thomas Peter^ Tim Romer^ 
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Prony’s method is a prototypical eigenvalue analysis based method for the re¬ 
construction of a finitely supported complex measure on the unit circle from its 
moments up to a certain degree. In this note, we give a generalization of this 
method to the multivariate case and prove simple conditions under which the 
problem admits a unique solution. Provided the order of the moments is bounded 
from below by the number of points on which the measure is supported as well 
as by a small constant divided by the separation distance of these points, stable 
reconstruction is guaranteed. In its simplest form, the reconstruction method con¬ 
sists of setting up a certain multilevel Toeplitz matrix of the moments, compute 
a basis of its kernel, and compute by some method of choice the set of common 
roots of the multivariate polynomials whose coefficients are given in the second 
step. All theoretical results are illustrated by numerical experiments. 
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1 Introduction 

In this paper we propose a generalization of de Prony’s classical method [10] for the parameter 
and coefficient reconstruction of univariate finitely supported complex measures to a finite 
number of variables. The method of de Prony lies at the core of seemingly different classes 
of problems in signal processing such as spectral estimation, search for an annihilating filter, 
deconvolution, spectral extrapolation, and moment problems. Thus we provide a new tool to 
analyze multivariate versions of a broad set of problems. 

To recall the machinery of the classical Prony method let C* := C \ {0} and let fj G C* 
and pairwise distinct G C*, j = 1,... ,M, be given. Let 6zj denote the Dirac measure in 
Zj on C* and let 

M 

i=i 
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be a finitely supported complex measure on C*. By the Prony method the fj and Zj are 
reconstructed from 2M + 1 moments 



Since the coefficients pg G C, i = 0,..., M, of the (not a priori known) so-called Prony 
polynomial 

M M 

p{Z) ■.= ^{Z-z,) = Y,PlZ^ 
j=i t=o 

fulhll the linear equations 

M M M M 

pef{£ -m) = Y m = 0,..., M, 

e=o j=i (.=0 j=i 

and are in fact the unique solution to this system with pM = 1, reconstruction of the Zj is 
possible by computing the kernel vector (pi,... ,pM-i, 1) of the rank-M Toeplitz matrix 

T:= (/(fe-^))^=o,...,M 

k=0^...,M 

and, knowing this to be the coefficient vector of p, compute the roots zj of p. Afterwards, the 
coefficients fj of / (that did not enter the discussion until now) can be uniquely recovered 
by solving a Vandermonde linear system of equations. When attempting to generalize this 
method to finitely supported complex measures on Cf, it seems natural to think that the 
unknown parameters Zj G Cf could be realized as roots of d-variate polynomials, and this 
is the approach we will follow here. As in the univariate case, the coefficients will be given 
as solution to a suitably constructed system of linear equations. However, for d > 2, an 
added difficulty lies in the fact that a non-constant polynomial always has uncountably many 
complex roots, so that a single polynomial cannot be sufficient to identify the parameters 
as its roots. A natural way to overcome this problem is to consider the common roots of a 
(finite) set of polynomials. These sets, commonly called algebraic varieties, are the subject 
of classical algebraic geometry and thus there is an immense body of algebraic literature on 
this topic from which we need only some basic notions as provided at the end of Section 2. 

Our main results are presented in Section 3, which is divided into three parts. In Sec¬ 
tion 3.1 we prove sufficient conditions to guarantee parameter reconstruction for multivariate 
exponential sums by constructing a set of multivariate polynomials such that their common 
roots are precisely the parameters. In Section 3.2 we focus on the case that the parameters 
Zj lie on the d-dimensional torus, which allows us to prove numerical stability, provide some 
implications on the parameter distribution, and construct a single trigonometric polynomial 
localized at the parameters. Finally, we state a prototypical algorithm of the multivariate 
Prony method. 

In Section 4 we discuss previous approaches towards the multivariate moment problem 
as can be found in [18, 1] for generic situations, in [27, 25, 17] based on projections of the 
measure, and in [7, 6, 4, 5] based on semidefinite optimization. Finally, numerical examples 
are presented in Section 5 and we conclude the paper with a summary in Section 6. 
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2 Preliminaries 


Throughout the paper, the letter d gN always denotes the dimension, C* 
let 

d := {Cd = {z G & : Zi ^ 0 foT alU = 1,... ,d} 


C \ {0}, and we 


be the domain for our parameters. For 2: G C^, k G Z'^, we use the multi-index notation 
:= z\^---z^‘^. We also letT ;= {z G C ; \z\ = 1} and the d-fold Cartesian product 
is called d-dimensional torus. We start by defining the object of our interest, that is, 
multivariate exponential sums, as a natural generalization of univariate exponential sums. 


Definition 2.1. A function /: Z*^ —)■ C is a d-variate exponential sum if there is a finitely 
supported complex measure p, on Cf, such that for all k G Z'^, f{k) is the k-th moment of p, 
that is, if there are M G N, fi,, fu G C*, and pairwise distinct zi,..., zm £ such that 
with p := fj^zj have 


f{k) 



M 

x^dp{x) = fjZ^j 


i=i 


for all k G 

In that case M, fj, and zj, j = are uniquely determined, and f is called M- 

sparse, the fj are called coefficients of f, and Zj are called parameters of f. The set of 
parameters of f is denoted by Qf or, if there is no danger of confusion, simply by fi. 

Remark 2.2. Let fj G C* and pairwise distinct tj G [0,1)'^, j = 1,... ,M, be given. Then 
the trigonometric moment sequence of t = Ylj=i 


f: Z'^^C, 


k 1-^ 



M 

i=i 


(where kt denotes the scalar product ofk and t) is a d-variate exponential sum with parameters 
g27ntj _ j'g27ritj,i^ _ _ ^ g yjpi analyzed in detail in Section 3.2. 

Let /; Z*^ —7- C be an M-sparse d-variate exponential sum with coefficients fj G C* and 
parameters Zj G C^, j = 1,...,M. Our objective is to reconstruct the coefficients and 
parameters of / given an upper bound n for M and a finite set of samples of / at a subset of 
Z*^ that depends only on n, see also [23]. 

The following notations will be used throughout the paper. For n G N, let In := {0,..., n}^ 
and let N := ]/„] = {n-\- 1)'^. The multilevel Toeplitz matrix 

Tn{f) := 

which we also refer to as T^, will play a crucial role in the multivariate Prony method. Note 
that the entries of are sampling values of / at a grid of \In — In\ = (2n -|- l)'^ points. 

Next we establish the crucial link between the matrix and the roots of multivariate 
polynomials. To this end, let 

n := C[Zi, ...,Zd] = {Y^PkZ’l^ ■■■Zf-.Fcni finite, pu G C}. 

k&F 
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denote the C-algebra of d-variate polynomials and for p = ^ n \ {0} let 

maxdeg(p) := max{||A:||oo : Pk / 0}. 

The A^-dimensional subvector space of d-variate polynomials of max-degree at most n is 
denoted by 

Iln := {p G n \ {0} ; maxdeg(p) < n} U {0} = spanjC'^ 3 z : k e In}, 
and the evaluation homomorphism at n = {zi,..., zm} will be denoted by 

, pr3 {p{zi),... ,p{zm)), 

or simply by An- Note that the representation matrix of An w.r.t. the canonical basis of 
and the monomial basis of n„ is given by the multivariate Vandermonde matrix 

A _ (r- fpAfxTV 

k^In 

The connection between the matrix and polynomials that vanish on Q lies in the observa¬ 
tion that, using Dehnition 2.1, the matrix Tn admits the factorization 

Tn = {f{k - i))k,eel„ = PnAlDnAn, (2.1) 

with Dn = diag(d), dj = j = 1,... ,M, and a permutation matrix Pn G {0, 

Therefore the kernel of An, corresponding to the polynomials in n„ that vanish on fl, is a 
subset of the kernel of T„. 

In order to deal with the multivariate polynomials encountered in this way we need some 
additional notation. The zero locus of a set P C 11 of polynomials is denoted by 

V{P) := {z G ; p{z) = 0 for all p G P}, 

that is, V{P) consists of the common roots of all the polynomials in P. For a set kl C C'^, 

7(17) ;= {p G n : p{z) = 0 for all z G 17} = hei A^ 

neN 

is the so-called vanishing ideal of 17. Finally, for a set P C 11 of polynomials, 

m 

{P) ■= {X] m G N, G U, pj G Pj 
i=i 

is the ideal generated by P. Note that V{P) = V{{P)) always holds. Subsequently, we identify 
Tin and C'^ and switch back and forth between the matrix-vector and polynomial notation. 
In particular, we do not necessarily distinguish between An and its representation matrix An, 
so that e.g. “V(ker makes sense. 

3 Main results 

In the following two subsections, we study conditions on the degree n, and thereby on the 
number (2n -|- 1)^ of samples, such that the parameters Zj can be uniquely recovered and the 
polynomials used to identify them can be computed in a numerically stable way. 
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3.1 Complex parameters, polynomials, and unique solution 

Our first result gives a simple but nonetheless sharp condition on the order of the moments 
such that the set of parameters 0 and the zero loci V(keT: An) and VikeiTn) are equal. 

Theorem 3.1. Let /; Z'^ —)■ C 6e an M-sparse d-variate exponential sum with parameters 
Zj G j = 1,..., M. If n> M then 

nf = V{kevTn{f)). 

Moreover, if this equality holds for all M-sparse d-variate exponential sums f, then n > M. 

Proof. Let Ll := Llf = {zi,... ,zm}- We start by proving Ll = y(ker^„). Since T(ker^„) D 
V{An+i) L) n, it is sufficient to prove the case n = M. It is a simple fact that I{{zj}) = 
{Zi — Zj^i,..., Zd — Zj^d), and that these ideals are pairwise comaximal, and hence we have 

M M 

= iYliZij - Zj,ej) Aj G {1,... ,4) C (ker^M) C /(O), 
i=i j=i 

which implies {ker Am) = I{Ll). Thus we have ^(ker^A^) = k^((ker^A^)) = ^ 

where the last equality holds because fl is finite (and can easily be derived from the above). 

Thus it remains to show that kei Am = kerTjvf. We proceed by proving rank^M = M. To 
simplify notation, we omit the subscript M on the matrices. Let N := dimllM and suppose 
that A G has rank r < M. Let Ll' = {zi,..., Zr} and w.l.o.g. let A' G denoting 

the first r rows of A, be of rank r. Now the first part of the proof implies the contradiction 
n = Vikei A) = Vikei A) = Q.'. 

Considering the factorization T = PAADA as in Equation (2.1) and applying Frobenius’ 
rank inequality (see e.g. [15, 0.4.5 (e)]) yields 

XcvakPA D + rank HA — rankH < rank HA = rankT < rank A 

which implies rankT = rank A = M. The factorization clearly implies ker A C kerT which 
together with the rank-nullity theorem dim ker A = N — M = dim kerT yields the hnal result. 

The converse follows from the fact that for 17 := {{xj, !,...,!) G Clf : j = 1,... ,M} with 
distinct Xj G C*, any subset B C n„ such that V{B) = 17 (which holds, by assumption, for 
B = kerT„) necessarily contains a polynomial of (max-)degree at least M. ■ 

Example 3.2. Let f be a 3-sparse 2-variate exponential sum with parameters Zj G and 
17 = {zi, Z 2 , zs}. The generating system of I{II) given in the proof of Theorem 3.1 is 

3 

P ■.= {pt ■. t e {1, 2}^} , Pii^i, Z 2 ) := Y\{Zej - Zj/.). 

i=i 

We start by illustrating the generic case that no two coordinates are equal, i.e., Zj/ A Zi/ if 
J / 1 and i = 1,2. The zero locus of each individual polynomial p^ is illustrated in Figure 3.1, 
where each axis represents C. The zero locus of each linear factor is a complex curve and 
illustrated by a single line. We note that the set P is redundant, i.e. the last three polynomials 
in the first row of Figure 3.1 are sufficient to recover the points uniquely as their common 
roots, but there is no obvious general rule which polynomials can be omitted. 
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Four other point configurations are shown in Figure 3.2. In the first three configurations 
coordinates of different points agree, which allows to remove some polynomials from P. In 
particular, the third point set which is collinear is generated already by {Zi — zij)(Zi — 
Z 2 ,i)(Zi — 2 : 3 , 1 ) and (Z 2 — 2 : 1 , 2 )^- The fourth point set is generated either by the above set P 
of polynomials or by Y^j^fiZi + Z2 — Zj^i — 2:^,2) and Z\ — Z2, (which are not elements of P). 


Figure 3.1: Zero sets of the polynomials in P, d = 2, M = 3, and for the case that no two 
coordinates are equal. 



Figure 3.2: Point sets Q C C^, d = 2, M = 3. 

Remark 3.3. Concerning the “natural” generator P = {n" liZij - Zj^g.) : G {1,..., d}| 
used in the proof above, we note that although the ideals {P) = (ker^iv^) coincide, the sub¬ 
vector space inclusion 

span P C ker Am 

is strict in general as can be seen for d = 2, M = 2, zi = (0,0), 2:2 = (1,1) and the 
polynomial Z\ — Z 2 G ker^ 2 - Moreover, we have the cardinality |P| = d^, at least for 
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different coordinates Zj^(, and thus |P| » — M = divahei Am, i-^-, the generator P 

contains many linear dependencies and is highly redundant for large M. 

Finally, we would like to comment on the degree n and the total number of samples fln+lY 
with respect to the number of parameters M: 

i) A small degree nGN, M<N<M + d, and surjective An results in an uncountably 
infinite zero locus V(kei An), since dim(ker^„) < N—M < d and thus is generated 
by less than d polynomials. 


a) Increasing the degree results “generically” in a finite zero locus, cf. [1], but “generically” 
identifies spurious parameters since e.g. for d = 2 Bezout’s theorem yields \ V{p,q)\ < 
deg(p) deg(g) with equality in the projective setting (counting the roots with multiplicity), 
for coprime polynomials p,q £ keiAn- 

Remark 3.4. We discuss a slight modification of our approach. Instead of In = {0,..., = 

{/c G Nq ; ||A:||j^ < n} we take Jn '.= {A: G Ng : ||A:||^ < n} as index set and consider the matrix 


Hnif) ■■= ifik + i))k,ieJ„ G 


instead ofTnif) = {f{k—i))k/£i„. Theorem 3.1 also holds with Tn replaced by Hn with almost 
no change to the proof. In this way we need only rather than (2n + 1)'^ samples of f 

and also allow for arbitrary parameters Zj G C'^ instead of Zj G Cf. While Tn is a multilevel 
Toeplitz matrix, Hn is a submatrix of a multilevel Hankel matrix, and for the trigonometric 
setting discussed in the following subsection, it is more natural to consider the moments f{k), 
k G ||A;||oo < n, than f{k), k G Ng, ||A:||^ < 2n. 


3.2 Parameters on the torus, trigonometric polynomials, and stable solution 

We now restrict our attention to parameters Zj G T*^, hence Zj = for a unique tj G 

[0,1)'^. In this case, VikevAn) fulfills a 2'^-fold symmetry in the following sense. Let p{z) = 
J2k=oPkZ^ £ ker^n and z = {zi,..., Zd)~^ £ with p{z) = 0, then P = (z\~^,Z 2 ,... ,Zd)^ 
is a root of the Ist-coordinate conjugate reciprocal polynomial 

n 

q{z) := pC"p(W-^,Z 2 , ...,Zd) = '^Pn-ki,k 2 ,-,kA^ ■ 

k=0 

Since the roots z £ Q. <Z are self reciprocal / = z, we have q £ ker^„ and thus z G 
y(ker^„) implies z' £ y(ker^„) for all choices of a conjugated reciprocal coordinate. 
Moreover, we have the following construction of a so-called dual certificate [7, 6 , 4, 5]. 

Theorem 3.5. Let d,n,M £ N, n > M, tj £ [0,1)*^, j = 1,...,M, zj := and 

Q := {zj : j = I,..., M} be given. Moreover, let p£ £ , i = 1,..., N, be an orthonormal 

basis with pe, £ ker(T„)-*-, I = 1,...,M, and pe : C, pe{z) = Pe,kZ^, then 

p : [0, ^ C, 

1 ^ 

?(*) = p-i) 

' e=i 

is a trigonometric polynomial of degree n and fulfills 0 < p{t) < 1 for all t £ [ 0 , 1 )*^ and 
p{t) = I if and only if t = tj for some j = 1,..., M. 
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Proof. First note that every orthonormal basis p£ G , £ = 1,..., A^, leads to 


N 




-2)1^ = 


e=i 


N N N 

^ Fz® '^pe,rPe,s = 

r,5=l 


IFP = N 


e=i 


r=l 


for z G T'^. Moreover, z = z ^ on T'^ yields that p is indeed a trigonometric polynomial. 
Finally, Theorem 3.1 assures J2f=M+i \p£{z)\‘^ = 0 if and only if z G fl. ■ 

We proceed with an estimate on the condition number of the preconditioned matrix T = Tn- 

Definition 3.6. Let M G N and Ll = ; tj G [0, 1)“^, j = 1, ..., M}, then 

sep(n) := min \\tj—t£ + r\\oo 

is the separation distance ofLl. For q > 0, we say that is (/-separated i/sep(fl) > q. 

Theorem 3.7. Let d,n,M G N, tj G [0, 1)'^, j = 1, , M , Zj := g > 0, and Q := {zj : 

j = 1,...,M} be q-separated. Moreover, let fj > 0, then n > 2dq~^ implies the condition 
number estimate 

mnd,WTW < . max, f, 

where the diagonal preconditioner W = diagtc, Wk > 0, k G In, is well chosen. In particular, 
lim„_^oo cond2imF = maxj fj/ miuj fj. 

Proof. First note, that the matrix T is hermitian positive semidefinite and define the condition 
number by cond 2 r := ||T|| 2 11^112, where Tl denotes the Moore-Penrose pseudoinverse. Let 
D = diagd, dj = f^"^, J = 1, • • •, M, and K = AW^A* G , then we have 


cond2lTTVF = cond2W A* AW = coud^DAW^ A* D = 


2 A* 


max fj 
min fj 


cond2iL 


and Corollary 4.7 in [20] yields the condition number estimate. The second claim follows since 
lim„_s.oo cond2iF =1. ■ 


In summary, the condition 

n > max{2(i(/“^, M} (3.2) 

allows for unique reconstruction of the parameters and stability is guaranteed when com¬ 
puting the kernel polynomials from the given moments. 

Remark 3.8. Up to the constant 2d, the condition n > 2d/q in the assumption of Theorem 3.7 
is optimal in the sense that equidistant nodes tj = j/m, j G Im, n < q~^ = m, imply 
A G C™' and rankyl = = M. We expect that the constant 2d can be improved 

and indeed, a discrete variant of Ingham’s inequality [19], [27, Lemma 2.1] replaces 2d by 
C^fd but gives no explicit estimate on the condition number. 

Moreover, Theorem 3.1 asserts that the condition on the degree n with respect to the number 
ob parameters M is close to optimal in the specific setting. We briefly comment on the 
following typical scenarios for the point set and the relation (3.2).- 



i) quasi-uniform parameters tj G [0,1)^ might he defined via sep(17) ps CdM i.e., 
ma,x{2dq~^, M} = M, 

a) equidistant and eo-linear parameters, e.g. tj = imply sep(i7) Ri CM~^, 

i.e., both terms are of similar size, 

in) and finally parameters tj G [0, l)'^ ehosen at random from the uniform distribution, 
imply Esep(r2) = , see e.g. [29], and thus max.{2dq~^, M} = 

Dropping the condition n > M in (3.2) and restricting to the torus, we still get the following 
result on how much the roots of the polynomials in the kernel of T can deviate from the original 
set Q. 

Theorem 3.9. Let d,n,M G N, tj G [0,1)'^, j = 1,..., M, Zj := q> t), and D := {zj : 

j = 1,..., M} be q-separated, then n > 2dq~^ implies 

D C T(kerT) n T'^ C : z G D, ||t||oo < 2d/n}. 

Proof. We prove the assertion by contradiction. Let y G T(kerT) n be 2(i/n-separated 
from the point set D C and let pi G C'^, ^ = 1,... ,N — M, constitute a basis of ker T. By 
definition, we have pe{y) = 0 and thus the augmented Fourier Matrix 

■ 4 , ■=(()). 

fulhlls AyPi = 0, i.e., dim ker > N — M. On the other hand. Corollary 4.7 in [20] implies 
rankyly = M + 1 and thus the contradiction N = dim ker + rankylj^ > + 1. ■ 

3.3 Prototypical algorithm 

Let / be an M-sparse d-variate exponential sum with pairwise distinct parameters Zj G Cf 
and n > M be an upper bound. Theorem 3.1 justifies the following prototypical formulation 
of the multivariate Prony method. 


Algorithm 1 Multivariate Prony method. 
Input: d, n G N, 

f{k), k G {-n,...,nY 
Set upr„ = (/(A:-£)),,,g,„GC^x^ 

Compute kerT„ 

Compute V{kei Tn) 

Output: VfkeiTn) = {zi,... ,zm} 


The third step, i.e., the computation of the zero locus V (ker T„), is beyond the scope of this 
paper and several methods can be found elsewhere, see e.g. [2, 22, 33, 34]. We further note 
that the number (2n + l)'^ of used samples scales as O (M'^) and that standard algorithms 
for computing the kernel of the matrix have cubic complexity. 
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4 Other approaches 


There are many variants of the one dimensional moment problem from Section 1, originating 
from such diverse fields as for example signal processing, electro engineering, and quantum 
chemistry, with as widespread applications as spectroscopy, radar imaging, or super-resolved 
optical microscopy, see e.g. the survey paper [24]. Variants of Prony’s method with an in¬ 
creased stability or a direct computation of the parameters without the detour via polynomial 
coefficients include for example MUSIC [32], ESPRIT [30], the Matrix-Pencil method [16], the 
Approximate Prony method [26], the Annihilating Filter method [35], and methods relying 
on orthogonal polynomials [12]. 

Multivariate generalizations of these methods have been considered in [18, 1] by realizing 
the parameters as common roots of multivariate polynomials. In contrast to our approach, 
both of these papers have an emphasis on the generic situation where e.g. the zero locus 
of two bivariate polynomials is finite. In this case, the total number of used moments for 
reconstruction might indeed scale as the number of parameters but no guarantee is given for a 
specific instance of the moment problem. A second line of multivariate generalizations [27, 25] 
decomposes the multivariate moment problem into a series of univariate moment problems 
via projections of the measure. While again this approach typically works well, the necessary 
number of a-priori chosen projections for a signed measure scales as the number of parameters 
in the bivariate case [17]. We note that the subset 

M 

Po := lUiZe - Zj,e) d}c I{n), 

i=i 

of the set of generators in the proof of Theorem 3.1 are exactly the univariate polynomials 
when projecting onto the d coordinate axes, see also the first and last zero locus in Figure 3.1. 

A different approach to the moment problem from Section 1 has been considered in [7, 6, 
4, 5] and termed ‘super-resolution’. From a signal processing perspective, knowing the first 
moments is equivalent to sampling a low-pass version of the measure and restoring the high 
frequency information from these samples. With the notation of Remark 2.2 the measure r 
with parameters tj G [0, l)'^ is the unique minimizer of 

min||z/||TV s.t. [ diz{t) = f{k), k e Im 

J[o,ip 

provided the parameters fulfill a separation condition as in Section 3.2. This is proven via 
the existence of a so-called dual certificate [7, Appendix A] and becomes computationally 
attractive by recasting this dual problem as a semidefinite program. The program has just 
(n -|- 1)^/2 variables in the univariate case [7, Corollary 4.1], but at least we do not know 
an explicit bound on the number of variables in the multivariate case, see [11, Remark 4.17, 
Theorem 4.24, and Remark 4.26]. 

Finally note that there is a large body of literature on the related topic of reconstructing a 
multivariate sparse trigonometric polynomials from samples, see e.g. [3, 21, 13, 8, 28, 31, 14]. 
Translated to the situation at hand, all these methods heavily rely on the fact that the 
parameters tj G [0,1)'’* are located on a Cartesian grid with mesh sizes 1/mi,..., l/nid for 
some mi,..., rrid G N and deteriorate if this condition fails [9]. Hence, these methods lack one 
major advantage of Prony’s method, namely that the parameters tj G [0,1)'^ can, in principle, 
be reconstructed with infinite precision. 
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5 Numerical results 

All numerical experiments are realized in MATLAB 2014a on an Intel i7, 12GByte, 2.1GHz, 
Ubuntu 14.04. 



(a) Sum of squared absolute val- (b) Zero set of sum of squared ab- (c) Squared absolute values of 
ues of kernel polynomials on T solute values of kernel polynomi- the first polynomial orthogonal 
(identified with [0,1)). als on C. to the kernel. 


Figure 5.1: Parameters d = 1, M = 3, n = 30, ti = 0.12, t 2 = and Dashed 

lines and x indicate no weighting, solid lines and • indicate triangular weights 
Wk = min{A: + 1, n — k}, A: = 0,..., n — 1. 


Example 5.1 [d = 1). We consider the case d = 1 with parameters on the 1-torus T that 
we identify with the interval [0,1). For a 3-sparse exponential sum some of the associated 
(trigonometric) polynomials are visualized in Figure 5.1, where we start with the upper bound 
n = 30 > 3 and also indicate the effects of a preconditioner W according to Theorem 3. 7 on 
the roots of the polynomials. 

The method introduced in [7] finds a polynomial of the form (3.1) as a solution to a convex 
optimization problem, whereas we find such a polynomial with Prong’s method. For this com¬ 
parison we used the MATLAB code provided in [7] and modified it so that it runs for different 
problem sizes depending on the sparsity M = 1,..., 100. This means that we used roughly 
5M samples and random parameters tj G [0,1), j = 1,... ,M, satisfying the separation con¬ 
dition in [7]. We only measured the time for finding a polynomial of the form (3.1), since 
the calculation of the roots is basically the same in both algorithms. In Figure 5.2 (a), where 
the times needed with cvx are depicted as circles and the times needed by Prong’s method 
are depicted as crosses, we see that the solution via convex optimization takes considerably 
more time. Note that the end criterion of the convex optimization program is set to roughly 
10“®, therefore the solution accuracy does not increase beyond this point, whereas for Prong’s 
method the solutions in this test are all in the order of machine accuracy, 10 “^®. 


Example 5.2 {d = 2). We demonstrate our method to reconstruct the parameters from the 
moments /: —>■ C, /c i—)• (1,1)^ + (—1, —1)^. For moments of order \k\ < n = 2 and the 

associated space of polynomials II 2 with reverse lexicographical order on the terms, we get the 
9x9 block Toeplitz matrix T = T 2 with the Toeplitz blocks T',T" as follows: 

/r T" T'\ /2 0 2\ /O 2 0\ 

T = if{k - £))k,e^i, = \T" r T"\, T' = 0 2 0 , T" = 2 0 2 . 

\r' T" T'J \2 0 2 / \0 2 0 / 
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(a) Running time for Prony method (+) (b) Running time for Prony method (+) 
resp. super-resolution (o) for varying number of resp. super-resolution (o) for varying number of 
parameters in the case d = 1. parameters in the case d = 2. 

Figure 5.2: Time comparison for extracting a polynomial of form (3.1), once with cvx depicted 
as circles and via Prony’s method depicted as crosses. 


A vector space basis of her T is given by the polynomials 

Pl = -1 + Zf, P2 = -Zi + Z 2 , P3 = -1 + Z 1 Z 2 , 

Pa = ~Z\ + Z\Z2, P5 = —1 + Z 2 , Pq = —Zi + Z\Z\, 

Pi = + z\z\. 

Since p^ = pi + Zip 2 , Pa = Zips, p^ = Z 2 P 2 + Ps, P6 = Zip^, and pj = {1 + ZiZ2)p3, we 
have (kerT) = {pi,P 2 ) and hence P(kerT) = V{pi,p 2 ) = {(1,1), (—1, —1)}. The zero loci 
ofpi,p 2 are depicted in Figure 5.3 (a) (in the style of Figure 3.1) resp. (b), where the torus 
T2 is identified with [0,1)^. Note that we would typically expect the intersection of the zero 
locus of each polynomial with the torus to be finite, which is the case neither for pi nor p 2 . 
In Figure 5.3 (c) the sum of the squared absolute values of an orthonormal basis of kerT is 
drawn. 

Example 5.3 {d = 3). Figure 5.4 depicts the intersection o/T^ (identified with [0, and 
the zero loci of two polynomials that arise with the Prony method for M = 2 parameters 
choosing n = 1 (which is not an upper bound for M). This illustrates that, in the case d = 3, 
the zero locus of a single polynomial intersected with the torus can typically be visualized as 
a “one-dimensional” curve as suggested by the heuristic argument that a complex polynomial 
can be thought of as two real equations, which together with the three real equations that define 

qpS 

as a subset of provides five equations, thus leaving one real degree of freedom. 


6 Summary 

We suggested a multivariate generalization of Prony’s method and gave sharp conditions 
under which the problem admits a unique solution. Moreover, we provided a tight estimate 
on the condition number for computing the kernel of the involved Toeplitz matrix of moments. 
Numerical examples were presented for spatial dimensions d = 1,2,3 and showed in particular 


12 





if 

/ 


(a) See also Fig. 3.1, zero 
loci V{pi),V(j> 2 ) C 



0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 

(b) Zeros of the kernel polyno- (c) Sum of squared absolute val- 
mial Pi (lines) and p 2 (dashed) ues of kernel polynomials, 
on 


Figure 5.3: Parameters d = 2, M = 2, n = 2, ti = (0.0,0.0) and t 2 = (0.5,0.5). 



(a) Zeros of the kernel polynomial pi on T®. (b) Zeros of the kernel polynomials pi, P 2 on T^. 

Figure 5.4: Parameters d = 3, M = 2, n = 1, ti = (0.1,0.3,0.25) and t 2 = (0.7,0.8, 0.9). 
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that a so-called dual certificate in the semidefinite formulation of the moment problem can 
be computed much faster by solving an eigenvalue problem. 

Beyond the scope of this paper, future research needs to address the actual computation of 
the common roots of the kernel polynomials, the stable reconstruction from noisy moments, 
and reductions both in the number of used moments as well as in computation time. 
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